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ANALYTIC PREDICTION OF AIRPLANE EQUILIBRIUM 

SPIN CHARACTERISTICS 

By William M. Adams, Jr. 

Langley Research Center 

SUMMARY 

The nonlinear equations of motion are solved algebraically for conditions for which 
an airplane is in an equilibrium spin. Constrained minimization techniques are employed 
in obtaining the solution. Linear characteristics of the airplane about the equilibrium 
points are also presented and their significance in identifying the stability characteristics 
of the equilibrium points is discussed. Computer time requirements are small making 
the method appear potentially applicable in airplane design. 

Results are obtained for several configurations and are compared with other analytic - 
numerical methods employed in spin prediction. Correlation with experimental results is 
discussed for one configuration for which a rather extensive data base was available. A 
need is indicated for higher Reynolds number data taken under conditions which more 
accurately simulate a spin. 


INTRODUCTION 

Modern fighter airplanes tend to have inertia and aerodynamic characteristics 
which produce unsatisfactory handling qualities and stability near the stall that can lead 
to inadvertent entry into the post stall/spin region during maneuvers at high angles of 
attack. In some cases a steady developed spin can result unless correct recovery tech- 
niques are promptly initiated. More than 200 fighter and trainer airplanes were lost in 
post stall/spin accidents during the period from 1966 to 1970 at a cost of more than 
360 million dollars in airplanes and more than 100 fatalities (ref. 1). 

Experimental studies of the spin characteristics of airplanes involving tests of 
dynamically scaled models have been made at the Langley Research Center for a number 
of years and yield valuable information. The accuracy of the data obtained is limited, 
however, and some parameters of interest are not measured. Scale effects such as 
critical Reynolds number differences can make extrapolation of results to the full-scale 
airplane difficult. 



Analytical techniques offer the possibility of spin analyses early in the design 
process that could provide information complementary to that obtained experimentally. 
Analytical methods currently employed in spin prediction fall into two categories. One 
type gains the capability of at-the-desk determination of approximate solutions by either 
neglecting the effects of some parameters or requiring that only part of the conditions 
for an equilibrium spin be satisfied (refs. 2,3, and 4). The second type begins with con- 
ditions hopefully either near the spin or favoring development of a spin and integrates the 
equations of motion forward in time until the average values of the variables are approxi- 
mately constant (refs. 5 to 8). Such integrations can require considerable amounts of 
computer time, and residual oscillations in the variables usually necessitate estimation 
of the equilibrium conditions for steady spins. In some cases persisting oscillations 
seen with this method may allow identification of a spin as being oscillatory rather than 
steady. 

An analytical technique allowing algebraic determination of airplane equilibrium 
spin conditions is presented in this paper which is limited in accuracy only by the mass, 
inertia, and aerodynamic data available. Numerical integration is avoided by casting the 
problem in a form requiring a constrained minimization solution of a nonlinear function 
of several variables. Solutions are obtained by using a method for finding a local mini- 
mum developed by Davidon (ref. 9) and modified by Fletcher and Powell (ref. 10). Results 
are presented for several configurations and are compared with other analytical predic- 
tions. Correlation with experimental results is discussed for one configuration for which 
an extensive set of aerodynamic data was available. 

Linearized equations of motion are also developed to provide information about the 
stability of the equilibrium spin conditions. A limited investigation is made of the char- 
acteristics of the system about the equilibrium solutions. 

SYMBOLS 

Values are given in both SI and U.S. Customary Units. The measurements and cal- 
culations were made in U.S. Customary Units. 

A,B matrices in linearized state equation £ = A£ + Bu 

A v ,By,C v coefficients of a quadratic equation, defined in appendix B 

A X> A Y> A Z body axis components of aerodynamic forces divided by mV 

a x = Ax + G^, a^- = Ay + Gy, a2 = A2+G2 
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b 


wing span 


c mean aerodynamic chord 

F external force on airplane (vector sum of aerodynamic and 

gravitational forces) 

f vector of nonlinear functions yielding f = (o!,j3,p,q,r) 

G vector of nonlinear functions yielding = (df,4,V,p,q,r) 

Gy jGy ,G ^ body axis components of gravity force divided by mV 

g acceleration due to gravity (assumed constant), 

9.80665 m/sec 2 (32.174 ft/sec 2 ) 

h airplane altitude above earth's surface 

*X’*Y’*Z’*XZ body axis moments and product of inertia about the center of mass 

J function to be minimized in order to find an equilibrium spin condition, 

f T f 

Mx rolling moment acting about X body axis 

My pitching moment acting about Y body axis 

M 2 yawing moment acting about Z body axis 

m mass of airplane 

p,q,r angular rates about body axes 

1 2 

q^ dynamic pressure, -pV 

R radius of helical path of airplane 

A A A A. 

R,T,Zj unit vectors in cylindrical coordinates; Zj is directed toward 

center of earth (fig. 1) 
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s 


wing area 


T X’ T Y’ T Z 


t i>t 2 

t 

t* 

to 

U 

u 

u,v,w 

V 


torque components expressed in body axes and divided by the corresponding 
moment of inertia, TV = *-component of torque 

X ix 

matrices displaying coordinate system transformations (see appendix A) 
time 

time at which search for equilibrium solution is made 

initial or reference time in linearized equations of motion 

control vector, U T = (6 e ,6 a ,6 r j 

deviation from nominal control vector, U(t) - UnOO 

body axis components of V 

airspeed 


W weight 

X,Y,Z axes 

x,y,z position triple (when devoid of subscripts, a body-axis coordinate 

system is referred to) 

a angle of attack 

/3 angle of sideslip 

y flight -path angle, tan -1 ^, 

6 e> 6 a>^r elevator, aileron, and rudder deflections (positive 6 e is trailing edge 
down, positive 6 a is right trailing edge down, positive 5 r is 
trailing edge left) 


n T = (6,0,y,w,i//’) 
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Xi ,A.2>**»>^8 characteristic roots of linearized equations 

,e,4> 

v o 

£ deviation from nominal value of H 

p ' atmospheric density 

i^,6, <p angles defining transformation between inertial and body axes 

(see appendix A) 

\j/ angle from horizontal projection of airplane longitudinal axis to horizontal 

component of V (positive clockwise when looking up) 

u> angular rate about center of mass, \Jp% + q2 + r 2 

Subscripts: 

H horizontal component 

I referred to inertial axis system 

IND reading at boom attached to nose of model 

N evaluated on nominal flight path 

o evaluated at time t 0 

S stall angle 

w referred to wind axis system 

Special symbols: 

(") unit vector 

("*) vector 

( *) derivative with respect to time 
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transpose of matrix 


() T 


( f 


value which minimizes J 


O 


column vector 


Coefficients and derivatives: 

Results presented herein are given in terms of the aerodynamic coefficients and 
their derivatives defined in the tabulations that follow: 










Pitching moment 

Cm - 

My 


<loo SC 

Cm 6 e 

9 C m 

“ 96 e 

c m q 

_ 9 c m 

4 

8 qC 


2V 


Longitudinal force 

C X = 
Cx 6 e 

,F - mg) • x 

<loo S 

9C X 

36 e 


Z body axis 
force component 




PROBLEM FORMULATION AND SOLUTION 
Equilibrium Spin Requirements 

An algebraic method of solving for airplane equilibrium spin conditions which 
utilizes nonlinear programing techniques is discussed in this section. The solutions 
obtained satisfy the requirements of a steady developed spin if the equilibrium points are 
stable. Solutions that have linear representations that are unstable either correspond to 
oscillatory spin conditions or are not actual developed spin conditions. 

An airplane in an equilibrium spin has the following characteristics: 

(1) It is operating in the post stall region 

(2) Descent occurs in a helical path about a vertical axis 

(3) The axis of the helix, its radius, and the descent speed would remain constant 

were it not for atmospheric density variation with altitude and other 
disturbances 

Three coordinate systems are employed in this analysis: 

(1) An inertial system with origin on the axis of the helix at an altitude h 

(positive Z-axis downward) 

(2) A system fixed in the body with origin at the center of mass 

(3) A wind axis system having origin at the center of mass with X w alined 

with the airplane velocity vector 

The transformations connecting these coordinate systems are given in appendix A. 
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Having defined the necessary coordinate systems one can state the equilibrium con- 
ditions mathematically as follows. Figure 1 will aid in the interpretation of some of the 
terms. 



Figure 1.- Top view of helical path. 

o? > Q!g (Airplane is in post stall region) (1) 


cu = i//Zj 
V = Ri//T - hZj 




F = -mRoj 2 R 


(Helical motion persists) 


Expressing the forces in the wind axis system and solving for a, (3, and V in 
equation (5) allow equations (4) and (5) to be written in the form 
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Search for an equilibrium spin condition is made at a particular time point t* . 

The unconstrained airplane has six degrees of freedom with second -order equations of 
motion; therefore, 12 parameters are required to completely specify its state at a given 
instant. However each independent equality constraint imposed upon the system can 
remove one parameter. This fact is employed in the search for an equilibrium condition 
to reduce the number of parameters to be searched over. The search is made at a speci- 
fied altitude h* and is limited to regions where V = 0. In addition, search is made only 
in regions that satisfy the conditions stipulated in equations (2) and (3) by including these 
four equality constraints explicitly. Finally, the orientation of the horizontal inertial axes 
is assumed to be such that i//(t*) = 0. The dimension of the state space to be searched 
over is, therefore, reduced from 12 to 5 if each control surface is held in a fixed position. 

For convenience the set of variables to be searched over was chosen to be 
{0,(p,y , See figure 1 and the symbols for a definition of \jy . 

The equations of motion are presented in appendix B as are the consequences of 
explicitly satisfying the equality constraints at each point. Engine thrust terms are not 
included. Earlier studies have indicated that these effects tend to be small and incon- 
sistent in the spinning region; use of thrust in this region is normally avoided since seri- 
ous engine damage can result and flameouts are likely (refs. 11, 12, and 13). 

Solution for conditions satisfying the additional relations in equations (6) is required. 
Define 



From appendix B it can be seen that satisfying f = {V} requires the simultaneous solu- 
tion of five coupled nonlinear equations involving tabulated aerodynamic data. Only a 
numerical solution is possible and there is no guarantee that a solution will exist. The 
procedure employed in this report attempts to solve the following equivalent problem. 
Find rf such that 
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5 

t / \ V ,2 -2 o 2 -2 .2 .2 

J(t?) = > f i = a* + j3 +p +q +r 
i=l 

has a local minimum at t) = rf while requiring that search be limited to the subspace 
within which the following equality constraints are satisfied: 

\]4t*) = 0 

h = h* (A specified constant) 

u(rj) = i^Zj (i.e., oj X] . = w yj = 0) 

V(ri) = 0 /V(t?) is determined such that V(rj) = o\ 

\ (see appendix B) / 


If Of* > ag and J(t 7 *) = 0, an equilibrium solution has been found. 

As formulated for this problem, a solution is obtained in an iterative manner begin- 
ning with an initial estimate rjy The method of search is initially equivalent to a gradi- 
ent or steepest descent procedure, but second-order information is accumulated during 
the iterations and quadratic convergence is approached near the end of the iterations. 

Thus, the method has initially the sureness of convergence toward a local minimum of a 
gradient procedure. Terminally the rapidness of convergence of a second-order method 
is approached without the necessity of computing and inverting the matrix of second par- 
tial derivatives at each point. Details of the method are presented in references 9, 10, 
and 14. Reference 14 also includes an extension to inequality constraints and is the paper 
most nearly describing the methods employed in the computer program that was modified 
for this study. 

The procedure developed in this report will henceforth be called a function minimi- 
zation technique. The solution 77 * which minimizes J ( 77 ) is a local one. Additional 
equilibrium spin solutions may exist for the same control settings that would be found if 
a different initial estimate were chosen. In fact two solutions were usually found for a 
given control setting if one was found. These two solutions are identified as flat and 
steep spins with the flat spin having the more nearly horizontal longitudinal axis. Inverted 
spins were not investigated due to lack of aerodynamic data at negative angles of attack. 

The solutions found by using the function minimization technique involve no approxi- 
mations other than those present in the data. Therefore, given sufficiently realistic mass, 
aerodynamic, and moment of inertia data, the solutions obtained would represent actual 


V(77) = Ri//T - hZj 
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airplane equilibrium spin conditions. The stability of these solutions must be examined, 
however, to determine whether they represent developed spin conditions one might expect 
to approximately maintain in a model or flight test. 


System Linearized About an Equilibrium Spin Condition 

A linearized representation of the airplane equations of motion is developed in 
order to obtain some information about the stability of the equilibrium spin conditions. 
Density variation with altitude is neglected. Define E(t) such that 



a,P,rjr,VA,r,0,<t> 
v o 


where 

E(t) = G[s(t),U(tj] 

to first order 


i(t) = E(t) - * N (t) = (v e g)^ * + ju(t) - u N (t)] 


or 

|(t) = A(t) 4(t) + B(t) u(t) 


The elements of the matrices A and B are presented in appendix C. It is assumed 
subsequently that Ej^t) corresponds to an equilibrium spin condition with fixed nominal 
controls in which case A and B are constant matrices. 

All stable roots for the linear representation indicates a return to the equilibrium 
spin conditions for sufficiently small deviations therefrom. Additional analysis would 
be required to determine whether a solution with an unstable linear representation is an 
oscillatory or a no-spin solution. 


Mass, Inertia, and Dimensional Characteristics 

Sketches of the planforms of the configurations studied are shown as figure 2. Con- 
figuration A represents a swept -wing fighter, configuration B represents a delta-wing 
fighter, configuration C represents a supersonic trainer, configuration D represents a 
stub-wing research vehicle, and configuration E represents a twin -jet swept -wing fighter. 

The weight, inertia and dimensional characteristics of the various configurations 
are shown in table I. These data were taken from references 6, 7, 8, 13, and 15 to 21. 
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Configuration A 



Configuration D Configuration E 


Figure 2.- Plan views of configurations. 


TABLE I.- WEIGHT, INERTIA, AND GEOMETRIC CHARACTERISTICS 


Configuration 

Prediction 

method 

(a) 

W, 

N 

(lb) 

>x. 

kg -m2 

(slug -ft 2 ) 

l Y> 

kg-m2 

(slug-ft2) 

*Z> 

kg -m2 
(slug -ft 2 ) 

J xz* 

kg-m 2 

(slug -ft 2 ) 

S, 

m 2 

(ft 2 ) 

b, 

m 

(ft) 

c, 

m 

(ft) 

Center 
of mass, 
percent c 

Control deflection limits 


o e > de s 

6a> de E 

6 r > de S 


A 

S, I, E, M 

105 739 

15 875 

112 060 

120 985 

0 

35.80 

10.87 

3.61 

33 

-30 to 10 

±15 

±6 

7, 8, 15 



(23 771) 

(11 709) 

( 82 654) 

( 89 237) 

(0) 

(385.55) 

(35.67) 

(11.83 ) 






E 

F 

150 000 to 180 000 





50.01 

11.71 

4.8895 

















-21 to 9 

±30 

±30 




(33 000 to 40 000) 





(538.34) 

(38.41) 

(16.04 ) 






E 

S 

161 594 

35 397 

157 574 

178 457 

0 

50.01 

11.71 

4.8895 

33.9 

-21 to 9 

±30 

±30 

17 



(36 328) 

(26 108) 

(116 222) 

(131625) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 






E 

S 

156 221 

27 832 

152 225 

172 780 

0 

50.01 

11.71 

4.8895 
















29.6 

-21 to 9 

±30 


18 



(35 120) 

(20 568) 

(112 277) 

(127 438) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 






E 

S 

159 659 

30 674 

162 593 

184 478 

0 

50.01 

11.71 

4.8895 

32.4 

-21 to 9 

±30 

±30 

18 



(35 893) 

(22 624) 

(119 928) 

(136 066) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 






E 

S 

165 922 

28 148 

169 826 

193 146 

0 

50.01 

11.71 

4.8895 
















35.8 

-21 to 9 

±30 


18 



(37 301) 

(20 761) 

(125 259) 

(142 459) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 






E 

s 

154 064 

31 944 

161 994 

183 001 

0 

50.01 

11.71 

4.8895 


„1 t ^ 






(34 635) 

(23 561) 

(119 481) 

(134 975) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 






E 

D 













20 

E 

M 

160 968 

35 398 

157 576 

178 460 

0 

50.01 

11.71 

4.8895 








(36 187) 

(26 108) 

(116 222) 

(131 625) 

(0) 

(538.34) 

(38.41) 

(16.04 ) 

33.3 

-21 to 9 

±30 

±30 


B 

I, E 

110 365 

18 438 

173 539 

187 096 

0 

64.57 

11.62 

7.24 








(24 811) 

(13 600) 

(128 000) 

(138 000) 

(0) 

(695.05) 

(38.12) 

(23.755) 

30.0 

-25 to 10 

±7.5 

±25 

7, 8 

B 

M 

110 365 

18 438 

173 539 

187 096 

5884 

64.57 

11.62 

7.24 








(24 611) 

(13 600) 

(128 000) 

(138 000} 

(4340) 

(695.05) 

(38.12) 

(23.755) 

30.0 

-25 to 10 

±7.5 

±25 


C 

F 

56 270 





15.79 

7.70 

2.36 












— 




21.5 

-17 to 8 

±60 

±6 

13 



(12 650) 





(170 ) 

(25.25) 

( 7.73 ) 






C 

I, E, M 

44 653 

2 305 

39 995 

40 800 


15.79 

7.70 

2.36 








(10 038) 

( 1 700) 

( 29 500) 

( 30 100) 


(170 ) 

(25.25) 

( 7.73 ) 

21.5 

-15 to 5 

±60 . 

±6 

8 

D 

D, I, E, M 

55 936 

5 814 

99 462 

101 502 

0 

18.58 

6.82 

3.13 








(12 575) 

( 4 228) 

( 73 384) 

( 74 867) 

(0) 

(200 ) 

(22.36) 

(10.27 ) 

19.5 

-30 

±7.5 

±7.5 

6, 21 


F , flight test; D, drop test; S, spin tunnel; I, numerical integration; E, estimation method; M, function minimization. 


Aerodynamic data for configuration E were taken from references 22 and 23. Data 
representing the aerodynamics of configuration B were taken from reference 24. All 
other aerodynamic data were taken from reference 8 which contains a compilation of data 
for several configurations. Other references that present aerodynamic data for config- 
urations A, B, C, and D and discuss the manner in which some parameters were estimated 
include references 5, 6, 21, 25, and 26. 

The dynamic derivatives data used in this study contain several approximations: 

(1) The measured values taken for the dynamic derivatives include other inseparable 
terms ^e.g., C Zp actually is Cj p + sin a, Cn r is actually C nr - C n ^ cos a, etc.). 

See reference 23 for more discussion of the difficulties. 

(2) For all configurations other than E, constant values of C mn were used. These 

T. 

correspond to the values employed in reference 8 for configurations A, C, and D and to 
the value employed in reference 24 for configuration B. 


RESULTS AND DISCUSSION 


Equilibrium spin conditions predicted by using the function minimization technique 
are presented. These results are compared with other analytic predictions. Compari- 
sons are made between experimental test results and function minimization predictions 
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for configuration E and they indicate a need for additional aerodynamic data. A limited 
investigation is made of the stability of the equilibrium spin conditions predicted by func- 
tion minimization. 


Equilibrium Spin Predictions 

Experimental spin test results for configurations A, B, C, and D are presented in 
references 5, 6, 13, 15, 16, and 21. Explicit comparisons are not made herein between 
these experimental results and the function minimization predictions. The reader who 
makes such correlations is cautioned to give careful consideration to the fact that some 
of the aerodynamic data were estimated for these configurations possibly in such a man- 
ner as to favorably affect the correlations. 

In tables II to IV, spin condition predictions are compared which were obtained by 
using the function minimization technique, a numerical integration procedure, and the 
analytical estimation procedure presented in reference 4. The numerical integration 
procedure begins with conditions thought to be near the developed spin and then the equa- 
tions of motion are integrated forward in time until the average values of a and if/ 
are approximately constant. The estimation procedure assumes that spinning occurs at 
wings level flight (/3 = <p = 0) and requires only that q and r be zero and that the veloc 
ity be such that the drag balances the gravity force. 

Configuration A . - Table II indicates that, for configuration A, the average values of 
a, / 3 , and found by numerical integration correspond quite closely with those pre- 
dicted by function minimization. The smaller value of V predicted by numerical inte- 
gration is of little significance being due to a difference in atmospheric density. The 
integration was initiated at an altitude of 9144 m (30 000 ft), but equilibrium conditions 
were approached at some lower altitude where the atmosphere was more dense. Conse- 
quently a smaller velocity was required to obtain a balance between aerodynamic and 
gravity forces. 


TABLE II.- SPIN CHARACTERISTICS FOR CONFIGURATION A 
[Altitude, 9144 m (30 000 ft)] 


Prediction 

method 

6 e . deg 

6 a > de S 

6 r , deg 

a, deg 

P, deg 

7 . 

m/sec 

(ft/sec) 

co, rps 

— 
6, deg 

<t>, deg 

R, 

m 

(ft) 

deg 

Function 

-30 

5 

-6 

77.05 

-1.483 

84.70 

0.3900 

-12.92 

-0.8413 

0.3996 

-87.72 

minimization 






(277.9) 




(1.311 ) 


Numerical 

-30 

5 

-6 

73 to 81 

-6 to 3 

77.40 

.39 





integration a 






(254 ) 






Estimation 

-30 

5 

-6 

59 

0 

94.7 

.27 

-31 

0 









(311 ) 






Function 

-30 

15 

-6 

83.60 

-.9154 

82.94 

.5631 

-6.394 

-.6806 

.09778 

-87.77 

minimization 






(272.1) 




( .3208) 



a Reference 7. 
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The approximate analytic procedure predicted a lower spin rate and a much steeper 
spin with a higher value of V than that found by the function minimization technique. 

Configuration B.- Predicted spin characteristics of configuration B are shown in 
table III. Flat spin prediction by numerical integration beginning with an initial altitude 
of 9144 m (30 000 ft) is in good agreement with the function minimization prediction for 
a and /3 and yields a slightly higher value of u> (0.22 rps vs 0.1910 rps). The differ- 
ence in V between these two methods is again mainly due to the altitude loss that 
occurred during the time history generation. 

The analytic estimation method also predicts fiat spin characteristics that are near 
the function minimization values for a, to, and V. 

The steep spin predicted by the estimation procedure is slightly flatter and faster 
than that found by function minimization. Other steep spins predicted by function mini- 
mization are also shown. The spin rates are generally considerably lower than the flat 
spin rates. The zero control deflection solution shown is quite near the stall angle of 
attack and is highly unstable. 


TABLE HI. - SPIN CHARACTERISTICS FOR CONFIGURATION B 
[Altitude, 9144 m (30 000 It)] 


Prediction 

method 

Type ol 
spin 

6 e , deg 

«a. deg 

6 r , deg 

a, deg 

/3, deg 

V, 

m/sec 

(ft/sec) 


n 

0, deg 



Function 

Flat 

-25 

IB 

-25 

73.63 

-2.410 

78.36 

0.1910 

-16.40 


2.238 


minimization 







(257.1) 




( 7.342) 


Numerical 

Flat 

-25 


-25 

75 

-2.0 

66.8 

.22 





integration a 







(219 ) 






Function 

Flat 

-25 

1 

-25 

73.21 

-2.383 

78.55 

.1875 

-16.81 


2.366 


minimization 



1 




(257.7) 




( 7.763) 


Estimation 

Flat 

-25 


-25 

69 

0 

77.4 

.16 

-21 

0 










(254 ) 


• 




Function 

Flat 

0 

0 

0 

65.83 

-1.960 

79.40 

.1857 


.9331 

3.340 

-87.54 

minimization 







(260.5) 




(10.96 ) 


Function 

Steep 

-25 

■ 

-25 

49.56 

-2.978 

94.95 

.08234 

-38.92 

7.733 

28.93 


minimization 







(311.5) 




(94.98 ) 


Estimation 

Steep 

-25 


-25 

55 

0 

87.2 

.095 

-35 

0 










(286 ) 






Function 

Steep 

-25 


0 

49.18 

-3.161 

95.25 


-39.33 

7.665 

29.57 

-77.39 

minimization 







(312.5) 




EBB 


Function 

Steep 

-25 


-25 

59.15 

-1.225 

88.45 

.1239 

-30.24 

4.228 


-81.57 

minimization 







(290.2) 




(31.89 ) 


Function 

Steep 

0 


0 

35.59 

-9.395 

98.12 

.1743 

-54.35 

-3.772 

11.22 

-93.09 

minimization 







(321.9) 




(36.82 ) 



a Reference 7. 
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Configuration C.- Since the same aerodynamic model and moments of inertia, and 
so forth, were assumed in each case for configuration C, the numerical integration and 
function minimization predictions shown in table IV(a) should have been close. (The 
speed V was not given.) A relatively large difference is seen in spin rate however. 
From examination of figure 12 of reference 8 it appears that the spin rate could have 
still been building up when recovery controls were applied. Thus, with the numerical 
integration procedure, one can be mistaken about the location of a steady spin if a momen- 
tary hesitation in ip occurs while the general trend is toward increasing ip. To avoid 
such errors it would be necessary, generally, to perform the integrations for a longer 
period to be certain that steady conditions had been reached. 

The estimation procedure predicts a steeper spin (73° vs 83.71°) and a smaller spin 
rate (0.32 rps vs 0.5442 rps) than that found by function minimization. 


TABLE IV.- FLAT SPIN CHARACTERISTICS FOR CONFIGURATIONS C AND D 


Prediction 

method 

6 e , deg 

6 a > deg 

6 r , deg 

a, deg 

P, deg 

m/sec 

(ft/sec) 

of, rps 

8, deg 

<t>, deg 

R, 

m 

(ft) 

\p', deg 

(a) Configuration C; altitude, 12 192 m (40 000 ft) 

Function 

-15 

40 

-4 

83.71 

-1.116 

84.76 

0.5442 

-6.312 

-0.8103 

0.1343 

-94.17 

minimization 






(277.3) 




(0.4406) 


Numerical 

-15 

40 

-4 

81.0 

-5 to 1 


>.46 





integration a 












Estimation 

-15 

40 

-4 

73.0 

0 

87.2 
(286 ) 

.32 

-17 

0 



(b) Configuration D; altitude, 4572 m (15 000 ft) 

Function 

-30 

7.5 

-7.5 

88.0 

-2.339 

52.2 

0.4888 

-2.093 

-1.913 

0.1298 

-103.0 

minimization 






(171.2) 




(0.4260) 


Numerical 

-30 

7.5 

-7.5 

84 to 89 

1 to -8 

55.2 

.44 





integration ° 






(181 ) 






Estimation 

-30 

7.5 

-7.5 

85.0 


52.4 
(172 ) 

.32 

-5 

0 




a Reference 8. 
^Reference 7. 


Configuration D .- Flat spin predictions for configuration D are shown in table IV(b). 
Good agreement is seen between numerical integration and function minimization predic- 
tions of a, (3, and V and about a 10 -percent difference in spin rate. 

The estimation method predicts a much smaller spin rate than that found by function 
minimization (0.32 rps vs 0.4888 rps). 

Comparison of analytic methods . - The approximate method of reference 4 is fast 
and easy to use and requires only a desk calculator. It is probably more accurate for 
flat spins than steep spins. Because only part of the spin conditions are satisfied, the 
accuracy of a given solution would be difficult to ascertain without results from another 
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technique. Use of solutions obtained with this method as initial conditions for the numeri- 
cal integration procedure should reduce the time required for convergence to developed 
spin conditions. 

Spin prediction by numerical integration has several attractive features. Informa- 
tion is obtained about spin entry characteristics, and oscillatory spins can perhaps be 
distinguished from steady spins by using this method. Unless a good estimate of the 
spinning condition is available, however, considerable computer time may be required to 
converge upon the spin conditions. Also, a point where momentary hesitation in the time 
variation of spin rate occurs can be mistaken for a developed spin solution. 

The function minimization technique’s convergence to an equilibrium spin solution 
is relatively insensitive to the initial estimate. The data setup requirements are essen- 
tially the same as for a numerical integration. This is the only one of the analytic methods 
which precisely determines equilibrium conditions which allows constant coefficient linear 
equations to be written. Steady spin conditions are identified if the characteristic roots of 
the linear representation are stable. Additional study, possibly employing numerical 
integration, would be required to determine whether linearly unstable equilibrium spin 
conditions were oscillatory spins or not actual spins. However, beginning the numerical 
integrations with conditions only slightly perturbed from the equilibrium conditions should 
minimize the time required to make this determination. The computer time requirements 
for obtaining function minimization equilibrium spin solutions are small (on the order of 
2 seconds or less of Control Data 6600 central processor time per solution for multiple 
solution runs). This small time requirement makes the method appear potentially appli- 
cable as a design tool for use in evaluating the effects of aerodynamic parameter varia- 
tions on the spin characteristics of an airplane. 

Configuration E .- A number of experimental studies have been made of the spin 
characteristics of configuration E as can be seen in table V. Results are shown from a 
flight test and three spin tunnel model tests. These results are compared with function 
minimization predictions. Radio-controlled model drop tests are described in 
reference 20. 

Six sets of dynamic derivative data were available for configuration E, each of which 
was measured at one of two amplitudes of oscillation and one of three frequencies of 
oscillation (ref. 23). A set was chosen that was taken with the largest amplitude of oscil- 
lation and a reduced frequency of 0.156. With this rather arbitrary selection, function 
minimization predictions in table V are seen to be in close agreement with spin tunnel 
results for the flat spin. A similar agreement was observed between the drop test results 
and function minimization predictions made under conditions corresponding to those of the 
drop test. 
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TABLE V.- SPIN CHARACTERISTICS FOR CONFIGURATION E 

[Altitude, 7620 m (25 000 ft) except in flight tests where, e.g., in one of the flat spins 11 000 ra (37 000 ft) S h ^ 3048 m (10 000 ft)j 


Prediction 

method 

Type of spin 

6 e> 

deg 

6 a» 

deg 

»r, 

deg 

Center 
of mass, 
percent c 

a, 

deg 

0, 

deg 

V, 

m/sec 

(ft/sec) 

CO, 

rps 

e, 

deg 

deg 

R, 

m 

(ft) 

deg 

Reference 

Flight test 

Flat 

-21 

Variable 
(0 to 15) 

-30 

29 to 33 

82 

5 4 ^IND < 25 


0.22 





16 

Spin tunnel 

Flat 

-21 

30 

-30 

30.4 

84 


83.8 

.54 


-8 to 10 



17 









(275 ) 







Spin tunnel 

Flat 

-21 

30 

-30 

33.9 

85 


81.4 

.62 


1 



17 









(267 ) 







Spin tunnel 

Flat 

-21 

30 

-30 

29.6 

86 


79.6 

.54 


-4 to 6 



18 









(261 ) 







Spin tunnel 

Flat 

-21 

30 

-30 

32.4 

87 


79.6 

.63 


-6 to 8 



18 









(261 ) 







Spin tunnel 

Flat 

-21 

30 

-30 

35.8 

88 


79.6 

.60 


0 



18 









(261 ) 







Spin tunnel 

Flat 

-21 

30 

-30 

31.9 

83 


79.9 

.54 


+4 



19 









(262 ) 







Function 

Flat 

-21 

30 

-30 

33.3 

81.58 

-0.1541 

89.18 

.5244 

-8.694 

0.2351 

0.2227 

-125.1 


minimization 








(292.6) 




( .7307) 



Function 

Flat 

0 

0 

0 

33.3 

82.47 

.05369 

87.81 

.5903 

-7.804 

.3185 

.1562 

-136.1 


minimization 








(288.1) 




.5126) 



Spin tunnel 

Flat 

-21 

0 

-30 

30.4 

80 


87.2 

.47 


-6 to 5 



17 









(286 ) 







Spin tunnel 

Flat 

-21 

0 

-30 

33.9 

85 


87.2 

.50 


-4 to 3 



17 









(286 ) 







Spin tunnel 

Flat 

-21 

0 

-30 

31.9 

82 


79.4 

.52 


-6 to 6 



19 









(262 ) 







Function 

Flat 

-21 

0 

-30 

33.3 

82.24 

.1450 

88.75 

.5469 

-8.030 

.5000 

.1997 

-127.6 


minimization 








(291.7) 




( .6552) 



Flight test 

Steep 

=-21 

Variable 

=-30 

29 to 33 

60 



.11 to .17 





16 

Spin tunnel 

Steep 

-21 

30 

-30 

31.9 

46 to 58 


104 

.19 


-12 to 14 



19 









(342 ) 







Spin tunnel 

Steep 

-21 

30 

-30 

29.6 

58 to 71 


89.9 

.19 to .25 


-8 to 10 



18 









(295 ) 







Function 

Steep, non- 

-21 

30 

-30 

33.3 

47.99 

.9764 

109.4 

.2273 

-41.77 

5.939 

4.609 

-84.59 


minimization 

equilibrium 







(358.8) 




(15.12 ) 



Function 

Steep 

0 

0 

0 

33.3 

59.39 

-.1513 

93.42 

.2725 

-30.63 

2.288 

2.019 

-89.98 


minimization 








(306.5) 




( 6.623-) 



Function 

Intermediate 

-21 

30 

-30 

33.3 

73.17 

.05453 

90.98 

.3648 

-17.01 

1.049 

.6690 

-100.5 


minimization 








(298.5) 




( 2.195 ) 



Function 

Intermediate 

0 

0 

0 

33.3 

73.18 

-.4839 

89.28 

.3925 

-17.04 

.3416 

.5316 

-105.5 


minimization 








(292.9) 




( 1.744 ) 




The steep spin conditions predicted by function minimization lie within the range of 
oscillations noted in the experimental tests. However, the steep spin solution for control 
deflections of elevator up, rudder with, and aileron against the spin is a nonequilibrium 
one ^i.e., StJ (77) y/c)77 but J (77*) # oj. A time history beginning with the non- 

equilibrium conditions resulted in an oscillatory spin that tended to approach conditions 
resembling a flat spin. 

A spin mode intermediate to the flat and steep spins was also predicted by the func- 
tion minimization procedure and is shown in table V for two control settings. Such a spin 
mode for this configuration was noted in the spin tunnel tests described in reference 18. 

Subsequent runs with parameters from the other sets of dynamic derivative data 
generally did not correlate as well with the model test results. The results obtained are 
shown in table VI. This table indicates that the small amplitude, low angular rate, dynamic 
derivative data produced by forced oscillation of the model about static positions do not 
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TABLE VI.- DEPENDENCE OF SPIN CHARACTERISTICS UPON FREQUENCY AND AMPLITUDE 
OF OSCILLATIONS EMPLOYED IN DYNAMIC DERIVATIVE TESTS 


Reduced 

frequency 

Amplitude of 
oscillation, 
deg 

<*> 

deg 

0, 

deg 

m/sec 

(ft/sec) 

rps 

8, 

deg 

<t>, 

deg 

R, 

m 

(ft) 

deg 

0.109 

A</> = ±5 At// = ±5.25 

86.77 

0.9885 

87.17 

0.9002 

-3.496 

1.077 

0.07526 

-161.9 





(286.0) 




(0.2469) 


.156 

A</> = ±5 At// = ±5.25 


No flat equilibrium spin found to the right 


.203 

A <#> = ±5 At// = ±5.25 

83.58 

.2801 

87.66 

.6403 

-6.678 

.4943 

.1279 

-141.1 





(287.6) 




( .4196) 


.109 

A </> = At// = ±10.50 

82.80 

.2614 

87.72 

.6042 

-7.472 

.5120 

.1493 

-138.2 





(287.8) 




( .4899) 


.156 

A tp = At// = ±10.50 

82.47 

.05369 

87.81 

.5903 

-7.804 

.3185 

.1562 

-136.1 





(288.1) 




( .5126) 


.203 

A <p = At// = ±10.50 

79.91 

-.3203 

88.06 

.5090 

-10.38 

.07232 

.2322 

-125.9 





(288.9) 




( .7619) 



adequately represent the effects of the large amplitude, high angular rate spinning motion 
upon the aerodynamically produced forces and moments. This view is supported by ref- 
erence 17 which presents limited experimental rotary -balance test results showing that 
the aerodynamic forces and moments can be very nonlinear functions of the rate of rota- 
tion. Clearly, means of obtaining data which properly exhibit these nonlinear relation- 
ships are needed. Rotary -balance testing, described in reference 27, is designed to mea- 
sure aerodynamic forces and moments with the model rotating at a constant rate. These 
measurements can be taken with a nonzero spin radius and are made at fixed values of 
a and 0. Data taken by this technique might allow more reliable analytic prediction of 
model spin characteristics. 

Spin angles of attack predicted by model tests correspond closely with those observed 
in flight tests as do those predicted by function minimization, with the arbitrary set of 
dynamic derivative data discussed earlier being used. The same correspondence does not 
hold in spin rate however. The spin rates observed in the flight tests were much less than 
those found by model testing and function minimization for both flat and steep spins. Some 
differences are to be expected due to the limited accuracy of model and full-scale test 
procedures and the inadequacies in the aerodynamic data discussed previously. The large 
differences seen in table V, however, indicate the possibility that other error sources, 
such as a critical Reynolds number difference between model and full-scale conditions, 
are present. The differences are significant since the spin rate is one of the important 
factors determining the success or failure of spin recovery techniques. Reference 28 
describes techniques that have been employed in the Langley spin tunnel to partially com- 
pensate for effects of Reynolds number differences. Aerodynamic data obtained in high 
Reynolds number facilities are essential for analytic purposes if force and moment data 
exhibit significant Reynolds number dependence in the range between model and full-scale 
tests. 
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TABLE VIL- CHARACTERISTIC ROOTS OF LINEAR REPRESENTATIONS 


Configuration 

h, 

m 

(ft) 

6 e , 

deg 

a a> 

deg 

deg 

Type of spin 

a, 

deg 

( 0 , 

rad/sec 

Aj, *2> sec 

A 3 , A 4 , sec - * 

* 5 . * 6 , 6ec ’ 1 

\rj, sec - * 

Ag, sec'* 

A 

9 144 
(30 000) 

-30 

5 

-6 

Flat 

77.05 

2.454 

(-0.128, ±3.26) 

(-0.0119^2.76) 

(-0.110, ±2.44) 

(-0.206,0) 

(-0.0545,0) 

B 

9 144 
(30 000) 

0 

0 

0 

Flat 

65.83 

1.167 

(-0.169, ±2. 57) 

(-0.0476, ±1.62) 

(-0.121, ±1.16) 

(-0.147,0.0364) 

(-0.147,-0.0364) 

B 

9 144 
(30 000) 

0 

0 

0 

Steep 

35.59 

1.095 

(1.34 ,±1.35) 

(-0.959 ,±1.63) 

(-0.0957 ,±1.10) 

(-2.51,0) 

(-0.198,0) 

B 

9 144 
(30 000) 

-25 

7 

-25 

Flat 

73.21 

1.178 

(-0.161, ±2.73) 

(-0.0665, ±1.75) 

(-0.122, ±1.17) 

(-0.218,0) 

(-0.0705,0) 

B 

9 144 
(30 000) 

-25 

7 

-25 

Steep 

49.56 

.5174 

(-0.350 ,±2.35) 

(-0.141, ±0.908) 

(-0.103, ±0.518) 

(-0.177,0) 

(2.39 x 10“ 7 ,0) 

C 

12 192 
(40 000) 

-15 

40 

-4 

Flat 

83.71 

3.419 

(0.183, ±3.84) 

(-0.0687, ±3.42) 

(-0.296, ±3.38) 

(-0.226,0) 

(-0.0349,0) 

E 

7 620 
(25 000 ) 

0 

0 

0 

Flat 

82.47 

3.709 

(-0.0232, ±4.02) 

(-0.0895, ±3.68) 

(-0.484, ±3.41) 

(0.365,0) 

(-0.223,0) 

E 

7 620 
(25 000) 

0 

0 

0 

Intermediate 

73.18 

2.466 

(0.332, ±2.94) 

(-0.245, ±2.89) 

(-0.111, ±2.47) 

(-0.631,0) 

(-0.229,0) 

E 

7 620 
(25 000) 

0 

0 

0 

Steep 

59.39 

1.712 

(0.2 13, ±2. 48) 

(-0.626, ±2. 10) 

(-0.101, ±1.71) 

(-0.182,0) 

(0.155,0) 

E 

7 620 
(25 000) 

-21 

30 

-30 

Flat 

81.58 

3.295 

(-0.123, ±3.85) 

(-0.0996,±3.28) 

(-0.433, ±3.22) 

(0.464,0) 

(-0.219,0) 

E 

7 620 
(25 000) 

-21 

30 

-30 

Intermediate 

73.17 

2.292 

(-0.124,±2.77) 

(0.1 58 ,±2.74) 

(-0.109, ±2.29) 

(-0.543,0) 

(-0.227,0) 

E 

7 620 
(25 000) 

-21 

30 

-30 

Steep, non- 
equilibrium 

48.00 

1.428 

(0.405,±2.38) 

(-0.808 ,±1.82) 

(-0.0897 ,±1.43) 

(-0.140,0) 

(3.06 x 10'®, 0) 


Linear Characteristics 

A general development of the linearized equations of motion is presented in appen- 
dix C. In this section the linearized representation is utilized to obtain some indication 
of the stability of the equilibrium spin conditions. A complete stability analysis, which 
would involve consideration of the stability of nonlinear systems, has not been attempted. 

Characteristic roots for the system linearized about equilibrium spin conditions 
are shown in table VII for several of the configurations studied. A return to the equilib- 
rium conditions is to be expected for sufficiently small initial disturbances from stable 
equilibrium points. Figure 3 illustrates this type of behavior for the flat spin of config- 
uration B corresponding to control settings of elevator up, aileron against, and rudder with 
the spin. The time history is generated by using the nonlinear equations of motion and 
by assuming constant atmospheric density and an initial perturbation of Aq Q = -7.5°/sec. 
The oscillatory behavior and the damping characteristics seen in figure 3 correspond 
closely, particularly for a, to those of the second largest complex conjugate pair of 
characteristic roots for this equilibrium spin solution as given in table VII. This indi- 
cates that the mode associated with this pair of roots is the dominant one under a pertur- 
bation in pitch rate. 

Most linear representations shown exhibit one or more unstable roots as seen in 
table VII and require additional analysis to determine whether the equilibrium point cor- 
responds to a persisting oscillatory spin or whether it is not a developed spin condition. 
Such an analysis has not been attempted herein; however, limited results are shown in 
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Figure 3.- Time history for configuration B. = -7-5°/sec. 

figure 4 for the flat spin solution found for configuration E for zero control deflections. 
These time histories are generated by assuming constant atmospheric density and initial 
conditions perturbed from the equilibrium point by Aq Q = -5°/sec. 

The solid curves in figure 4(a) are generated by using the nonlinear equations of 
motion whereas the dashed curves were obtained by using a linear representation. Two 
sets of linear responses are shown for each parameter, one beginning on the nonlinear 
trajectory at t = 0 with the linear representation about the flat spin and the second 
beginning on the nonlinear trajectory at t = 10 seconds with the linear representation 
about the intermediate spin. Good qualitative agreement is seen for approximately 
10 seconds in each case indicating that the linearized equations of motion may yield an 
adequate qualitative description of the system's time development over a relatively large 
region around an equilibrium point. 

The results shown in figure 3 and figure 4(a) indicate that the linearized stability 
characteristics of the equilibrium spin conditions can yield valuable insight into the non- 
linear motions of an airplane in the vicinity of these points if the correct aerodynamic 
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Nonlinear equations of motion 

Linearized equations of motion 



(a) £q t=0 = -5%ec. 

Figure 4.- Time histories for configuration E. 


data in the spin regime can be obtained. Unfortunately the data are often inadequate. 

For example, flight and model tests of configuration E have demonstrated that the flat 
spin is stable and virtually steady. This result is in sharp contrast with the flat spin 
found by using the mathematical model as is seen in figure 4(a) and table VII. 

Figure 4(b) illustrates the nonlinear system time development with the same initial 
conditions as in figure 4(a) for C nr = -0.177 when a > 80°. The motion, although not 
steady, diverges much more slowly from the flat spin equilibrium point and illustrates 
the sensitivity of the dynamic characteristics to the dynamic derivatives. The value of 
C nr chosen corresponds approximately both to its value at the equilibrium point and to 
the rotation balance equilibrium value at a = 85° as determined from reference 17. 

The dynamic behavior of an airplane is strongly dependent upon its linear stability 
characteristics. However, other factors must also be considered if a desired dynamic 
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Time, sec 

(b) Ai 0 = -5%ec; C n ^ = -0.177 for 80 ° $ a $ 90 °. 

Figure 4.- Concluded. 

response is to be obtained. Thus a linearly unstable equilibrium point does not neces- 
sarily indicate that a recovery can be readily achieved. Some factors that could delay or 
prevent recovery include (1) other equilibrium points may exist in the vicinity, (2) the 
recovery controls applied may not excite the unstable mode or modes sufficiently, and 
(3) nonlinear effects may cause the spin to be an oscillatory one which tends to remain 
within a bounded region in the vicinity of the equilibrium point. One example of this type 
of behavior is a recovery attempt from the zero control flat spin of configuration E with 
the control settings often referred to as "optimum recovery controls" (i.e. , elevator up, 
aileron with, and rudder against the spin). Recovery was not achieved after 10 turns 
despite the tendency to diverge from the equilibrium conditions as discussed previously. 
The same recovery procedure for the zero control flat spin of configuration B, which is 
stable, resulted in a recovery in less than two turns; this recovery could be anticipated 
since no equilibrium spin condition was found for this configuration for full aileron with 
the spin. 

These examples indicate that caution should be exercised when attempting to cor- 
relate the stability characteristics of the linearized representation with ease of recovery. 


23 



Nevertheless, more detailed analysis of the factors affecting the stability of the equilib- 
rium spin conditions might aid in the design of configurations having more acceptable spin 
characteristics and identify more destabilizing recovery techniques allowing improved 
recovery capability. Reference 29 discusses some of the parameters affecting the sta- 
bility of very flat spins. 


CONCLUDING REMARKS 

A function minimization method has been demonstrated which will determine exact 
equilibrium spin conditions for any given inertia and aerodynamic data in the spin region. 
Linear characteristics of the equilibrium spin conditions are also obtained. Limited 
analysis indicates that the stability and dynamic characteristics of the linearized repre- 
sentation can yield valuable insight into the nonlinear motion of an airplane in the vicinity 
of the equilibrium points if the correct aerodynamic data can be obtained. 

Application of the method to the prediction of spin characteristics for one configura- 
tion for which an extensive data base was available yielded results which correlated well 
with model tests only when an arbitrary choice was made between available sets of 
dynamic derivative data. This arbitrariness coupled with incorrect stability characteris- 
tics for this configuration's flat spin mode indicates a need for more accurate aerody- 
namic data taken under conditions which more nearly represent the spinning models. 

The spin rates observed in flight tests of this configuration did not agree with those found 
in model tests or with those using data derived from models. The differences are due in 
part to the limited accuracy of model and full-scale test procedures and to failure to 
simulate adequately the environment of the spinning model when obtaining aerodynamic 
data; the large differences seen, however, indicate the possibility that other error sources 
such as a critical Reynolds number difference between model and full-scale conditions 
are present. 

The method can be employed to aid in the prediction of spin characteristics of future 
airplanes in the design stage if appropriate analytically or experimentally determined 
aerodynamic data are available. Computer time requirements per equilibrium spin solu- 
tion are small in multiple solution cases making the method appear applicable in studying 
the effect of aerodynamic parameter variations upon the spin characteristics of an 
airplane. 

Comparison of the spin characteristics predicted by function minimization with those 
observed in flight and model testing can aid in determining the adequacy of the aerody- 
namic data in the spin region. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 15, 1972. 
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APPENDIX A 


COORDINATE SYSTEM TRANSFORMATIONS 

Three sets of coordinate systems are employed in the analysis. The transforma- 
tions relating the coordinate systems are presented as follows. 

(1) Inertial to body axes (angular orientation): 



where 0 (yaw) is counterclockwise rotation about inertial Zj-axis, 6 (pitch) is coun- 
terclockwise rotation about resulting Y-axis, and 0 (roll) is counterclockwise rotation 
about an axis parallel to X-axis. Figure 5 shows, in three stages, the angular rotations 



Yaw Pitch Roll 


Figure 5.- Inertial to body axes transformation. 

defining the transformation. Each view is from a point on the positive axis about which 
the rotation is taken toward the origin of coordinates. 

The matrix T ^(0,0,0) has the following form 

cos 0 cos 0 sin 0 cos 6 -sin 0 

cos 0 sin 9 sin 0 cos 1 p cos 0 cos 9 sin 0 

Tj = - sin 0 cos 0 + sin 0 sin 9 sin 0 

cos 0 sin 9 cos 0 sin 0 sin 9 cos 0 cos 9 cos 0 

+ sin 0 sin 0 - cos 0 sin 0 

25 



(2) Body to wind axes (see fi 



where a (angle of attack) and /3 
in figure 6 and 

cos 0 sin 0 

Tg = -sin 0 cos (3 

0 0 
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(sideslip) are wind incidence angles positive as shown 



APPENDIX B 


EQUATIONS OF MOTION AND COMPLEMENTARY RELATIONSHIPS 

The six -degree -of -freedom nonlinear equations of motion are expressed in the fol- 
lowing form: 


a = q + 


'q S e 
— — Cv - — sin 0 + r sin 8) sin a 
mV x V 1 


fa S g 

+ C z + — cos 8 cos 0 - p sin pj cos a 


sec p 


|3= - 


C Y - — sin 0) sin p + r 
mV x V ' 


q ^ „ 

cos a + I — — Cv + — cos 8 sin 0 ] cos j3 
' mV * V 1 


s 

C 7 + — cos 6 cos 0 I sin p - p 
mV ^ V 1 


sm a 


V _ Cv - — sin 0 1 cos a cos p + [ Qy + — cos 0 sin 0) sin p 

V \mV A V / \mV 1 V / 


/q S 

+ ( C Z+ | cos 0 cos 0] sin a cos 0 


P = 


lY_Iz IxZ "L + (l-?X 7 ^)^ Z pq + . 


h 


l x l z 


l z / r x 


q °° Sb /r + Ixz P ^ 

vr V n y 


c + Iz - Ix nr + Ix -Z-fr2 p 2 ) 

" Iv m I V P Iv 1 P ; 


* - q °° Sc 

!y 111 Iy ‘ : Y 


r = 


l XZ 2 l Y ' 1 


X ’ 


^xiz I Z 


pq + 




1_ 

D 


where D = 1 - 


x XZ 

J X J Z 
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The aerodynamic coefficients are written in a form consistent with the data avail- 
able for the most extensively tested model. Note that this results in a linear dependence 
upon control deflections. 

C X = C X K/3,U=0) + C x (ct,0) 6 e 

°e 

Cy = C Y (a,/3,U=0) + C Yg (a,/3)6 e + c y 6 (« ,/3)5 a 

U G 3 . 

+ Cy 6r (a,/3)6 r + ^:JCy p (a!)p + Cy r (a!)r] 

C Z = C Z («AN + c Z fi ( a ,P) 6 e 

°e 

C t = C^(q!,/ 3,U=0) + C Zg (a,0)S e + C t (a,® 6 a 
+ C Zgr (a,^)6 r + JL|Cjp(a!)p + Cj r (at)rJ 

Cm = C m (o!,^,U=0) + C m ^(a,0)5 e + ~r C m q(a)q 

C n = C n (a 5 /3,U=0) + C n5 (af,/3)6 e + C ng (a,/3)6 a 
+ C n6r («,/3)6 r + — jC np (Q!)p + Cn r (Q!)r| 

The general relationships between the Euler angle rates and the body axis angular 
rates are 




r— 


— i 


r *\ 

* 


0 

sin 0 sec 8 

cos 0 sec 8 


P 

] 

>= 

0 

COS 0 

-sin 0 

< 

q) 

$ 

V. J 


1 

tan 8 sin 0 

tan 9 cos 0 


r 

L J 


The requirement that the angular velocity be vertical leads to the following 
relationships: 

p = -0 sin 0 

q = xf/ cos 6 sin c p 
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r = if/ cos 9 cos cp 
I * t 

with c o = h// . 

Using the stipulated condition that i//(t*) = 0 and expressing the linear velocity in 
terms of the body axes as follows allow determination of a and p in terms of 9, (p, 
y, and \p': 


r -\ 




r 

V cos y cos 

u 


V cos a cos p 

1 

1 

>=< 

V sin p 

} = T 1 (^o,e,0)< 

V cos y sin if/ \ 

w 


V sin a cos p 


-V Sin y 

V > 


L J 

1 

v. J 


Explicit inclusion of V = 0 constraint is accomplished by solving for the speed at 
which V = 0 given the other variables. 

The resulting expression for V is 

-By + sgn A v \j By - 4AyCy 
' 2A^ 

where 


A v 


1 £S 

2 m 



cos a; cos P + sin a cos p + 


[c Y Kj3,U=0) 


+ C Y («,/3)5 a + C Y (<M)6 r 
°a °r J 



Bv = | |[c Y p(a)p + Cy r (a)rj sin p 

and 

C v = g[(cos 6 cos (p sin a - sin 9 cos a)cos p + cos 9 sin <p sin / 3 J 

If additional data should be taken which result in the need for a more general func- 
tional form for the aerodynamic coefficients, an analytic solution for the velocity at which 

V = 0 might not be possible. In this case it would probably be inadvisable to require that 

V = 0 at each point in the search. One would then increase the dimensions of the space 
to be searched over from 5 to 6 and include the term (V/V)^ in the function J that is 
to be minimized. 

The final equality constraints that are required to be satisfied at each point are 

V • R = 0 and Vh = R&>. This is accomplished by requiring that R = (V cos y)/w. 
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LINEARIZED EQUATIONS OF MOTION 

The composition of the matrices in the following linearized system of equations 
will be exhibited: 

i(t) = A(t)!(t) + B(t)u(t) 

where 

4 T (t) = f Aa ,A/3, ,Aq,Ar ,A0 

The following relationships allow the matrix elements to be written more compactly: 


<lcO 

-X 

q S 
q = —2— 
mV 

D = 

1 l xz 
r x x z 

V 

_ l Y ~ l Z 

ix 

T l Z " X X 

I Y 

J 3 = 

l X - Iy 

l z 

A x 

= qC x 

Ay = qCy 

A Z : 

- qc z 

G x 

= - — sin 8 
V 

Gy = ^ cos 8 sin 0 

G z : 

= ^ COS 9 COS (j) 

a x 

= A X + G X 

ay = Ay + Gy 

a z = 

: A z + G z 


qjsb 

q_.sc 


q^sb 

T x 

ii 

* ! 
O 

rp °0 /-» 

1 Y t '■'m 

X Y 

T z : 

x OO 

T n 

l Z 


where the coefficients are defined in appendix B. 

Consider one additional point concerning the notation which will avoid a prolifera- 
tion of subscripts in what is presented subsequently. It is to be understood by the reader 
that all partial derivatives written are to be evaluated on the nominal flight path. For 

example, is to be evaluated by employing values of the variables taken at time t 

on the nominal flight path. If the nominal flight path corresponds to an equilibrium spin 
condition, the partial derivatives will not, therefore, change with time. 
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First, the elements of matrix A are defined as follows: 


da 

da 


+ r sin 0^ + 


3Ar 


da 


cos a 


(a z - p sin 0) 


3A- 


X 


da 


sin a) sec 0 


3<i • ( 3Av 3A Z \ 

— = la - q)tan 0 - p cos a - r sin a + sin a + — — cos msec 0 

3/3 7 V 3/3 30 


(V/V Q } = ["( A X ' Gx) sin a + ( A Z~ G z )cos ajsec /3 


Ayy 

— = -cos a; tan 0 

3p 


da 

3q 


= 1 


dot 

— = -sin a tan 0 
dr 


da g 


30 V cos /3 

3o! __ -g 
dtp V COS 0 


(cos 8 sin a - sin 8 cos p cos a) 
cos 8 sin p cos a 


p cos a + r sin a + 
da 


1®X 


3A Z \ / 3A X \ 

sin a - cos alar, + 

3(3! / 


9a 


sin |3 + 


3A^ 

da 


3/3 (”/ ~ 1 dAy 

d^ = ~ ( a X C0S “ + a z Sin a ) cos 0 + a Y sin /3 j + COS 0 


1 3 A x 3A z \ 

— — — cos a + — — sin a sin 0 


30 


’( V Ao) 


(^■X “ Gx^cos a + ^A z - G z jsin a sin 0 


(Ay - Gy) - q ^r(C Yp P + Cy r r) 


COS 0 1 


cos 0 
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|U Sin or + q A c Yp cos f! 


»e = o 

aq 


If = - cos a + q ^ c Y r cos 0 


= ^jjsin 6 cos <p sin a + cos 6 cos a)sin 0 - sin 6 sin (p cos |Sj 


= i. c ? - g -f (sin (p sin a sin /3 + cos cp cos 8) 
o<p V 


9 vAo _ 

da 

W/Vp 


9A zVc 


9A 


X 


av sin a + 1 a„ + — — | cos a 

yX 8a J \ 2 da 


9A T 

„„ ■ = -la Y cos a + sin a - 
9/3 \ X Z 9/3 


9Ay 

cos /3 + - sin 8 

da 


dA v dAr 


sin /3 + a v + cos a + sin a cos 8 

\ Y 9)3 9/3 1 


dV 

dV 


= 2||a x cos a + A z sin a)cos /3 + Ay sin /3j - q ^( c Y p P + c Y r r ) sin 0 


9v/ v o - b ^ 

— — — = q — Cy sin /3 
9p 2V P 


9V/V 0 


8q 


= 0 


9V/V 0 - b ^ 

— (■ — - = q — Cy sin 8 
9r 2 V * r 


9 V V o an % 1 

— — — = - ^|J sin 0 cos <p sin a + cos 6 cos a)cos ft + sin 6 sin cp sin /3 

dv/v 0 


d <p V 


= - — cos 0(sin <p sin a cos /3 - cos <p sin j3) 


dg_ _ ( 97 x Jxz 97 z\ 1 _ 

9a \ 9a I x 9cn y D 
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_9p _ / ^X l XZ ^ Z ) \ 
W [dp I x dp ) D 


Op = 2 (t + ^ T \ _b_ 

8 (V/V 0 )"\ \ X Ix Z )~ I X 2V 


(V + C 'r r ) + if (<V + C "r r ) 


9 p 

9 P 


/, ^ T \ X XZ n A q oo Sb b ( r l XZ r 
(1 + J 3 ) q+ _to + — c„ t 


1 _ 

D 


iE 

9q 


Jr ^ r + ( 1+j3) fe P 

V 1 Ixiz/ 1 


1 _ 

D 


.& = 
9r 


, J 1 - 


! XZ L . q °° Sb b 




q + 


c, 


L xz 


I Y 2VT ir+ lv ^ ni 


1 _ 

D 


iE=o 

30 

±=0 

30 

8 £L = i!!l 

30 ! 30 ! 


3q _ 9T Y 

dp dp 


W V o) 


= 2 ( T Y ‘ | C ™q q J 


3q 


XZ 


= j r - 2 — -p 


3p 


lY 




q^sc ? 


3q I Y 2V m q 


Ji = j.p + 2^r 

3r z Iv 
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^=0 

80 


S£ = Ax Z ^X 

ao! I I z 3a 3a / D 

ar = / X XZ ^X + f£z\ 1_ 

30 \i z a/3 3/3 y D 


3r 

9 ( v /v 0 ) 



qopSb b 

I 7 2V 


2cz 

lx 


W P + C; 


3r 

3p 



q + 


q«,sb 

b / 

^xz 

*z 

2V\ 

Ax 



n 


+ Cr 


1_ 

D 


3r 

3q 


3r 

3r 



I 

ID 



1_ 

D 


3r 
3 6 


= 0 



a£ 

da 


= 0 



8g_ 

3V 


= 0 


3£ 

ap 


0 


80 

3q 


cos 0 


80 

ar 


-sin 0 


— =0 — = -q sin 0 - r cos 0 

30 00 H r 


(Cn p P + C nr r 
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i£=o 

^=0 

30 

3a 

dp 

3V 

1 

■2^ = tan 0 sin (p 

30 

3P 

3q r 

3r 


= q sec^0 sin <p + r sec^0 cos <p 

ae 


= tan 6 ( q cos 0 
30 v r 


Next, the elements of the matrix B are defined as follows: 


da 

36~ 


3a 

367 


= q^-Cx 6 sin a + C Zg cos ajsec p 


= 0 


- 22 . = 0 
36 r 


3/3 _ 
36 e q 


if” (CXfig cos a + ^Zg sin °^ sin 0 + Cy g COS j3" 


P ] 


= qC Yga cos p 
= qCy g cos /3 
*2>. 


^Cx g cos a + C Zg sin a^cos p + 


Cv. sin 


'] 


J - 

— = qCy* sin p 
°a 


= qCy sin P 
°r 

-fvfcs)* 


W D 
1 


r sin 0) 
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96 t 


*L« 

d6 a 


A = 0 
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96 t 
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( Ci6 r 1 
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Cm 6 e 
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